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Abstract 


The  contraction  mapping  method  for  frequency  estimation  in  the  presence  of 
noise,  identifies  the  cosine  of  the  frequency  to  be  detected  as  a  fixed  point  of 
a  certain  correlation  mapping.  At  its  heart,  the  method  provides  a  plan  for 
automatic  self  tuning  of  parametric  filters.  A  variant  of  the  method,  called 
the  HK  algorithm,  produces  recursive  zero-crossing  rates  (normalized  HOC 
sequences)  that  converge  to  the  frequency  of  interest.  A  statistical  explana¬ 
tion  for  the  contraction  mapping  method  as  epitomized  by  the  HK  algorithm 
is  provided  when  the  HOC  sequences  are  produced  by  bandpass  filters.  The 
outright  consistency  of  the  zero-crossing  rate  is  not  required.  Examples  show 
that  the  method  performs  quite  remarkably. 


Abbreviated  Title:  “Contraction  Mappings” 

Key  words  and  phrases:  Stationary,  Gaussian,  bandwidth,  fixed  point, 
parametric  filter,  convergence. 

AMS  subject  classification:  Primary  62A/10,  secondary  62M07. 


1  Introduction 


The  iterative  method  for  frequency  detection  discussed  in  He  and  Kedem 
(1989)  (HK),  Yakowitz  (1990)  (Y),  and  Kedem  and  Yakowitz  (1990)  (KY), 
locates  discrete  frequencies  by  iterations  of  a  certain  contraction  mapping.  If 

is  a  frequency  to  be  detected,  then  cos(cji)  is  a  fixed  point  of  the  mapping, 
and  the  method  produces  a  recursive  sequence  that  converges  to  the  fixed 
point.  Thus,  the  essence  of  the  method  is  the  location  of  fixed  points  by 
recursive  sequences.  In  general,  we  refer  to  this  method  as  the  contraction 
mapping  (CM)  method  .  A  specific  special  case  is  the  algorithm  suggested 
in  HK.  It  makes  use  of  the  zero-crossing  rate  observed  in  filtered  time  series. 

Suppose  {Cr,r  £  (—1, 1)}  is  a  parametric  family  of  time  invariant  linear 
filters  defined  from  a  real  valued  impulse  response  sequence  h(n;r),  n  — 
0,  ±1,  ■  •  which  satisfies  a  certain  requirement  about  r.  The  HK  algorithm 
is  given  in  terms  of  the  recursion 

r*+ 1  =  cos(7rJ 

where  7r  is  the  observed  asymptotic  zero-crossing  rate  after  the  application 
of  { Cr  }■  7 r  is  an  example  of  observed  (normalized)  higher  order  crossings 
(HOC).  If  u> i  is  the  frequency  to  be  detected,  then  under  appropriate  con¬ 
ditions,  r*.  converges  to  cos^x),  so  that  ^Tk  converges  to  W],  as  k  — ►  oo. 
Looked  at  from  a  more  general  point  of  view,  the  HK  algorithm  provides  a 
sequence  of  self-tuned  filters  whose  successive  application  yields  a  sequence  of 
convergent  zero-crossing  rates. 

In  this  paper  we  discuss  the  strong  consistency  of  the  HK  algorithm, 
and  illustrate  its  use  through  some  examples.  In  particular,  in  Example  6, 
the  algorithm  is  applied  in  the  search  of  a  diurnal  cycle  in  a  well  known 
precipitation  data  set.  Our  analysis  supports  the  hypothesis  that  tropical 
rainfall  contains  a  diurnal  cycle. 

After  a  short  introduction  to  the  CM  method  in  section  2,  we  illustrate 
the  method  in  section  3  by  appealing  to  two  useful  families  of  bandpass  filters. 
Section  4  discusses  the  statistical  basis  for  the  CM  method  as  epitomized  by 
the  HK  algorithm,  when  a  family  of  bandpass  filters  is  used  in  generating 
convergent  (normalized)  HOC  sequences.  The  main  result  is  Theorem  4.1, 
formulated  for  the  bandpass  case.  Basically  it  says  that  if  the  bandwidth  is 
made  to  shrink  at  a  certain  rate,  7rfc  — >  u-i  with  probability  one. 
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Our  approach  does  not  require  an  outright  consistency  of  the  zero-crossing 
rate.  Successive  narrowing  of  the  bandwidth  renders  the  zero-crossing  rate 
more  and  more  accurate. 

Experimental  results  in  section  4,  using  a  particularly  convenient  family 
of  bandpass  filters,  confirm  what  is  suggested  by  the  statistical  analysis, 
namely,  the  convergence  of  the  observed  (normalized)  HOC  sequences  toward 
the  desired  frequency  is  rather  fast  and  precise  as  the  bandpass  is  squeezed  in 
a  controlled  manner.  The  examples  show  that  the  CM  method  can  provide 
more  precise  estimates  than  does  periodogram  analysis. 


1.1  Definition  of  HOC 

Let  {Ztj,  t  =  0,  ±1,  ±2,  •  •  •,  be  a  real  valued  zero-mean  stationary  time  series, 
and  let  {Eg}  ,  6  E  0,  be  a  parametric  family  of  time  invariant  linear  filters. 
Denote  by  {Zt(9)}  the  filtered  series, 


Zt{9)  =  Ce{Z)t 

Let  J[, 4]  be  the  indicator  of  the  event  A ,  and  suppose  {Zt{9)}  is  a  real  valued 
process.  The  quantity 

N 

De  =  ]T%«(0)z«- i(fl)<o] 

t= 2 

gives  the  number  of  zero-crossings  (in  discrete  time)  observed  in 


Zi(9),  Z2(0),  •  •  • ,  Zn(9) 


The  family  {  Dg },  8  E  0  is  called  a  (observed)  higher  order  crossings  family,  or 
HOC.  The  corresponding  expected  HOC  family  is  given  by  {E(Dg)},0  E  0. 
The  sample  zero-crossing  rate  is  defined  from  the  normalized  HOC, 


70 


^ xDg 
N  -  1 


(1) 


When  {ZtJ  is  strictly  stationary,  the  asymptotic  zero-crossing  rate  is  defined 
by  the  almost  sure  limit  of  the  sample  zero-crossing  rate, 


7 e  =  lim  %  =  lim 

N—*oo  N—<-oo 


leDg 
N  -  1 


(2) 


Note  that  E(jg)  =  E( jg) 


2 


2  The  CM  method 


We  shall  adhere  to  the  mixed  spectrum  model 


Zt  =  A\  cos(u>i/)  +  B\  sin(uqt)  +  (t 


(3) 


where,  t  =  0,  ±1,  ±2,  •  •  •,  A\  and  B\  are  uncorrelated  normal  random  vari 
ables  with  E(Ai)  =  E(B\)  =  0,  and  Var(A\)  =  V ar(B\)  —  a2.  Further, 
suppose  {C<}  is  colored  stationary  Gaussian  noise  with  mean  0  and  variance 
<t|,  independent  of  the  Ai  and  B\.  The  noise  is  assumed  to  possess  an  abso¬ 
lutely  continuous  spectral  distribution  function  F^{u)  with  spectral  density 

/cM, u  e  [— «•,  «■]- 

Let  r  be  a  parameter  that  takes  values  in  (  —  1,1),  and  consider  the  para¬ 
metric  family  of  linear  filters  {£r}  defined  by  the  family  of  impulse  response 
sequences  {/i(n;r)},  n  =  0,  ±1,--*.  The  corresponding  family  of  trans¬ 
fer  functions  is  denoted  by  {H(u',r)},  u  €  [ — 7r,  7r].  Let  Zt(r)  and  ( fir ), 
t  =  0,±1,---,  be  the  filtered  process  and  filtered  noise,  respectively.  The 
essential  element  of  the  CM  method  is  the  assumption  that  the  real  part  of 
the  first-order  autocorrelation  of  the  filtered  noise  ( fir )  is  r.  In  symbols,  we 
assume  that  (HK,Y), 


r  =  3R 


EUr)Ct-i(r)]\ 

E\(fir)\2  ) 


where  the  overbar  denotes  complex  conjugate.  This  means  that  we  require 
the  fundamental  property , 


The  fundamental  property  (  4)  is  exhibited  by  many  filters  of  which  the 
AR(\)  filter  (also  known  as  exponential  smoothing) 

Ca  =  1  +  otB  -f  a2B2  +  •  •  •  ,  a  6  (— 1, 1)  (5) 


where  B  is  the  backward  shift,  is  the  clearest  example  when  used  in  conjunc¬ 
tion  with  white  noise  £t  (see  HK). 

Let  pfir)  be  the  real  part  of  the  first  order  autocorrelation  of  {Zt(r)}, 


Pi(r)  =  SR 


JS[Z«(r)5_i(r)n 

E\ Zt(r)\*  J 
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and  assume  that  Cr  “passes  u>i”.  That  is,  \H(u ;  r)|  is  positive  in  a  neighbor¬ 
hood  of  u\  for  some  r.  Then  (  4)  implies  that  (HK,  Y), 


Pi(r) 


f<712[|tf(-u1;r)|2  +  \H(ux,  r)|2]  cos(wj)  +  [f%  |%i;r)|Wc(w) 


+  |^;r)p]  +  \ H^-r^dF^u) 


(6) 


This  function  of  r  is  a  weighted  average  of  cos(u>i)  and  r.  It  admits  the  more 
compact  representation, 


p\(r)  =  r*  -f  C(r)(r  —  r*) 


where  r*  =  cos^),  and 


C(r) 


E\(t(r)\2 

E\Zt(r)\* 


(7) 

(8) 


Clearly  0  <  C(r)  <  1,  and  C(r)  =  1  when  the  filter  does  not  pass  u>\. 
When  C(r)  <  1  in  a  neighborhood  of  r*,  we  identify  the  mapping  (  7)  as  a 
contraction  at  r*  with  a  contraction  factor  C(r).  From  the  recursion 


nfe+i  =  Pi(rk) 


(9) 


we  obtain 


Pi(n)  =  r-  +  in  C(r,)](r  ~  O 

j=l 


If  n-=1  C(rj)  —>  0,  then 


rjt  -►  r 


(10) 


and 


cos  J(rfc)  — +  uii 


as  k  — >  oo.  We  now  summarize  the  preceding  discussion  in 


Theorem  2.1  Let  { Zt }  be  given  by  (  3),  and  let  { CT }  be  a  parametric  fam¬ 
ily  of  time  invariant  linear  filters  for  which  the  fundamental  property  (4) 
holds,  and  one  that  passes  u)\.  Let  p\{r)  be  the  real  part  of  the  first-order 
autocorrelation  of  the  filtered  process  Zt(r )  =  Cr(Z)t.  Fix  r0  £  (—1,1).  If 
rij=i  C(rj)  —>  0,  then  r defined  by  the  recursion  r^+i  =  px(rk),  converges 
to  cos(cui),  and 

cos -1(ri)  —4 


as  k 


oo. 
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The  recursion  (  9)  characterizes  the  CM  method.  It  suggests  the  proce¬ 
dure 

nt+i  =  pi(rk)  (11) 

where  pi(r)  is  a  suitable  estimator  of  p\(r)  from  a  time  series  Z\,  •  •  • ,  Zpj,  of 
size  N. 

In  the  real  Gaussian  case,  a  natural  estimator  is  obtained  from  the  “cosine 
formula”  (see  the  discussion  in  Barnett  and  Kedem  (1991)) 

Pi(r)  =  cos(£(7r)) 

and  the  algorithm  takes  the  specific  form  (HK) 

rk+ i=cos(7rJ  (12) 

where  7 Tk  is  the  observed  asymptotic  zero-crossing  rate  as  defined  in  (  2).  We 
can  now  see  that  the  HK  algorithm  (  12)  is  a  special  case  of  the  CM  method. 

In  the  Gaussian  case  there  is  a  certain  advantage  for  using  an  estimate 
from  zero-crossings.  In  the  mixed  spectrum  case,  the  sample  autocorrelation 
is  not  necessarily  a  consistent  estimator.  Likewise,  the  sample  zero-crossing 
rate  is  in  general  not  a  consistent  estimator  for  its  expectation  when  the 
spectrum  contains  jumps.  However,  the  asymptotic  zero-crossing  rate  falls 
with  probability  one  between  the  lowest  and  highest  nonnegative  frequencies 
in  the  spectral  support,  so  that  its  variability  can  be  controlled  by  the  ap¬ 
plication  of  bandpass  filters.  See  Kedem  and  Slud  (1991)  (KS)  for  a  detailed 
discussion  regarding  the  consistency  of  the  sample  autocorrelation  and  the 
sample  zero-crossing  rate. 

In  the  real  nonGaussian  case,  when  the  HOC  family  is  obtained  from 
bandpass  filters  with  narrow  bandwidth,  the  cosine  formula  provides  an  ex¬ 
cellent  approximation  to  the  first  order  autocorrelation,  and  the  HK  algo¬ 
rithm  (  12)  performs  remarkably  well.  Evidence  of  this  is  given  in  Examples 

2  and  6  below. 

3  Two  parametric  families  of  filters 

This  section  discusses  the  CM  plan  relative  to  two  useful  parametric  bandpass 
filters.  In  addition,  it  is  shown  that  the  contraction  factor  can  in  fact  be 
controlled  by  the  filter  bandwidth. 
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3.1  The  complex  exponential  filter 

A  useful  example  (Y,KY),  that  illustrates  the  contraction  mapping  algorithm 
(  9)  and  its  key  elements,  is  provided  in  terms  of  a  parametric  family  of  linear 
filters  defined  by  the  impulse  response 


h(n ;  r,  M )  = 


exp  ( in6(r ))/ \J2M  +  1 ,  \n\  <  M 
0,  jnj  >  M 


(13) 


where  6{r )  =  cos~1(2^'1r)  th  cos-1(r).  Because  of  the  particular  form  of 
this  impulse  response,  we  shall  refer  to  this  filter  as  the  “complex  exponential 
filter”.  The  corresponding  squared  gain  is 


\rr,  .,m2  1  sin  i  2M  +  I  t 

\H(w,r,M)\2  =  2M  +  i 


)) J 


%))] 


—  7 r  <  U)  <  7T  (14) 


which  we  recognize  as  the  Fejer  kernel  centered  at  0(r).  This  parametric 
family  also  motivates  a  certain  assumption  concerning  the  contraction  factor 
needed  for  the  proof  of  our  main  result  in  the  next  section. 

One  can  see  that  \H(u\  r,  M)\2  assumes  its  maximum  value,  which  is 
2 M  +  1,  at  u  =  6(r).  This  filter  is  essentially  a  bandpass  filter  centered  at 
0(r),  and  its  dominant  lobe  is  supported  over  the  interval 


J  =  (0(r)  -  2tt/(2 M  +  1),  0(r)  +  2ir/{2M  +  1)) 

For  sufficiently  large  M ,  we  can  think  of  /  as  the  effective  bandwidth  of  the 
filter.  The  complex  case  follows  the  footsteps  of  the  real  case  except  that 
now  we  use  the  real  part  of  the  concerned  autocorrelations.  Thus,  under  the 
assumption  of  white  noise ,  the  fundamental  property  (  4)  now  takes  the  form 

£[C,(r,  A/)]] 

£|C.(r,Af)P  j 

where  (t(r,M)  is  the  filtered  noise  (it  is  not  white  anymore).  It  is  easy  to 
check  that  (  15)  holds  for  |r|  <  2M/(2M  +  1).  Similarly,  the  correlation 
mapping  is  now  defined  by 

jE[Zt(r,M)Zt_1(r,M)}} 

Pl[)  X\  E\Zt(r ,  M)\2  j 
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where  Zt(r,M )  is  the  filtered  process,  and  the  dependence  of  pi(r)  on  M 
is  suppressed  for  notational  convenience.  Thus  defined,  pi(r)  =  pi(r,M) 
satisfies 


Pi(r)  -  r* 


E\Ct(r,M)\ 2 
E\Zt(r,  M)\2 


(r  -r")  =  C(r)(r-r‘) 


(16) 


which  is  a  contraction  mapping  at  r*  =  cos(u;i),  as  long  as  u>\  passes,  with 
contraction  factor  (in  the  present  complex  case  the  gain  is  not  symmetric; 
when  M  is  sufficiently  large  negative  frequencies  and  in  particular  —  uj\  are 
suppressed) 


2 

C{r)  =  - - ^ -  (17) 

where  again,  the  dependence  of  C(r)  on  M  is  suppressed  for  notational  conve¬ 
nience.  The  recursion  rk+1  =  pi(r^)  yields  the  convergence  rk  — >  r*,  provided 
nj=1C(rj)  -4  0,  as  k  — >  oo. 


3.1.1  Bounding  the  contraction  factor  by  the  bandwidth 

As  is  clear  from  the  preceding  paragraph,  enhanced  speed  of  convergence 
towards  r*  can  be  obtained  by  decreasing  the  contraction  factor  C(rj )  with 
each  iteration  of  the  CM  recursion  (KY).  This,  however,  can  be  achieved  by 
controlling  the  effective  bandwidth  of  the  filter.  The  following  result,  which 
helps  in  motivating  a  certain  assumption,  renders  this  idea  more  precise. 


Proposition  3.1  Under  the  assumption  of  white  noise,  let  r,  M)\2  be 

as  in  (  If),  and  assume, 

uq.  6  [wa,Wii]  C  X 


Define 


A  =  Ub  -  <X>a 


Then  there  exists  a  C(M),  which  does  not  depend  on  r,  satisfying 


0  <  C(r)  <  C(M)  <  1 


and  C(M )  is  of  order  0(A)  as  A  — >  0. 
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Proof:  By  centering  the  filter  at  0, 


r,  M)\2  =  —  0(r); 


2  M 


2M  +  1 


,M)  |2 


Therefore,  when  Ui,  0(r)  G  [u>a,u;i], 


I#  (A;  W 


Define 


2M  +  1 


C(M)  = 


2M  + 


-c2 


i2+^c2 


Then  in  light  of  (  17), 

0  <  C(r)  <  C{M)  <  1 
As  A  — >  0,  choose  M  =  0(1/A). 


(18) 


□ 


As  a  matter  of  fact,  by  minimizing  (  18)  with  respect  to  M,  the  optimal 
value  of  M  is  approximately  1.165/A. 

Proposition  2.1  can  be  demonstrated  graphically  as  well.  Figure  1  (a) 
shows  the  graph  of  pi(r,  20)  when  u>i  —  0.8  radians  per  unit  time,  and  a\  = 
=  1.  In  a  neighborhood  of  r*  =  cos(0.8)  =  0.697,  we  can  see  from  the 
figure  that  the  derivative  p\{r,  20)  is  very  close  to  0,  and  this  is  already  with 
M  —  20.  Therefore,  in  a  neighborhood  of  r*, 

0  «  /9j(r*, 20)  »  ^-(.r’20)-T.r*  =  C(r, 20) 
r  —  r* 

However,  the  effective  bandwidth  of  the  filter  in  a  neighborhood  of  0  =  0.8 
is  considerably  greater  than  0,  as  can  be  seen  from  Figure  1  (b). 


3.2  A  complex  filter 

He  and  Kedem  (1989)  considered  the  parametric  filter  defined  by, 

Zt(r;  M)  —  (l  +  eie^B)M  Zt 
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where  M  is  a  positive  integer,  r  G  (— 1, 1),  and  6{r)  G  (— x,7r).  For  sufficiently 
large  M ,  for  example  M  =  30,  this  filter  acts  as  a  bandpass  filter  with  impulse 
response  and  squared  gain  given  by,  respectively, 


and, 


h(n;  r,  M)  — 


M  \  e"M»,  n  =  0,  •  •  • ,  M 

n  I 

0,  otherwise 


\H(u-,r,M)\2  =4Mcos2M 


The  fundamental  property  is  obtained  by  defining, 


;:,c  os“i 

(- 

fill) 
2  ) 

|  cos(A)dF^(A) 

f\  cos™  | 

\dFc(\) 

(19) 


(20) 


(21) 


When  {C<}  is  white  noise  (HK), 

r  — >  cos(0(r)),  M  — ►  oo 


and  we  have  the  approximation  0(r)  «  cos-1(r).  As  in  the  previous  case,  it  is 
easy  to  see  that  this  approximation  holds  true  for  any  continuous  spectrum 
noise,  provided  the  spectral  density  is  sufficiently  smooth.  For  sufficiently 
large  M,  e.g.  M  =  100,  when  the  filter  is  centered  at  0(r)  near  u>\,  it  only 
passes  a  band  of  frequencies  in  a  neighborhood  of  (see  Figure  2  (b))  and 
we  obtain  from  (21)  the  mapping, 

pi(r)  =  cos(wi)  +  C(r)(r  -  cos(u;i)) 


with 

cm  = 

f|if(^;r,M)|*  +  £|C1(r)p 

In  Kedem  and  Lopes  (1991)  it  is  shown  that  there  exists  an  increasing  se¬ 
quence  { Mk }  such  that 

rk+  i  =  pi(rk,  Mk)  -*  cos^x),  k  -*  oo 
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Figure  2  (a)  shows  the  graph  of  pfir)  =  pfir,  100),  for  M  —  100,  u>  =  0.8, 
and  cr\  =  =  1.  As  with  the  complex  exponential  filter,  the  derivative  of 

pfir,  100)  near  cos(0.8)  is  close  to  0,  so  that  C(r )  =  C(r,  100)  is  rather  small, 
much  smaller  than  the  bandwidth  given  in  Figure  2  (b),  and  we  can  see  that 
pi(r)  is  a  contraction  at  r*  =  cos(u>i)  =  cos(0.8).  Furthermore,  the  figure 
gives  additional  evidence  that  the  contraction  factor  can  be  bounded  by  the 
bandwidth  in  a  neighborhood  of  r*. 

Remark  3.2.1  In  the  multiple  frequency  case,  it  was  proved  in  HK  that  in 
the  purely  discrete  spectrum  case  (er^  =  0),  for  fixed  0 ,  pfiO)  (obvious  nota¬ 
tion)  converges  to  the  frequency  closest  to  0,  as  M  — ►  oo.  However,  when 
white  noise  is  present,  the  convergence  is  toward  cos($).  In  other  words,  by 
shrinking  the  bandwidth  without  “sliding  the  filter” ,  the  frequency  closest  to 
the  center  of  the  frequency  is  lost.  The  frequency,  however,  is  not  lost  if  we 
decrease  the  bandwidth  and  slide  the  gain  toward  the  frequency  at  the  same 
time.  This  is  the  essence  of  the  CM  method. 


4  Convergence  in  the  bandpass  case 


Suppose  the  parametric  family  [CT]  defines  a  family  of  bandpass  filters  sat¬ 
isfying  (  4),  where  Cr  has  bandwidth  of  size  Wr  stretching  from  *ui  to  uj*  , 
[*<*;, a;*]  C  [0, 7r] .  Let  {Lb}  be  the  corresponding  HOC  family  obtained  from 
a  time  series  of  length  N.  Consider  the  asymptotic  zero-crossing  rate, 


7 r  =  lim  7r  =  lim  -r~r 
TV- oo  TV— oo  N  -  1 


From  (KS)  we  know  that  7r  is  most  likely  a  random  variable  (mixed  spectrum 
in  a  Gaussian  process),  but  that  it  (and  hence  also  its  expected  value)  falls  in 
the  bandpass  with  probability  one.  Thus,  by  controlling  Wr,  we  can  control 
the  variability  of  7r,  and  make  it  arbitrarily  small.  Because  this  fact  is  of 
vital  importance  in  the  following  development,  we  state  it  formally. 

Proposition  4.1  (KS).  Assume  that  the  real-valued,  stationary,  zero-mean 
Gaussian  process  { Ztj  has  normalized  spectral  measure  v  supported  on  J  U 
(— J)  where 

J  -  [*u;,u>*]  C  [0, 7r] 
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Then,  whether  v  has  atoms  or  not,  the  asymptotic  zero-crossing  rate  yr  lies 
in  J  with  probability  one. 

Thus,  7 T  falls  with  probability  one  between  the  lowest  (»o>)  and  highest  (a;*) 
nonnegative  frequencies.  As  a  consequence  also  E{ 7r)  falls  in  and 

we  have  the  following  fact. 

Corollary  4.1  Let  a  be  any  number  between  *)r  and  Efyr),  and  let  Wr  be 
the  length  ( bandwidth )  of  the  interval  Then 

l7r  -  E{~fr)\\  sin(a)|  <  Wr 


with  probability  one. 

The  purpose  of  the  corollary  is  to  motivate  Assumption  A2  below. 

The  Gaussian  assumption  and  the  contraction  mapping  (7)  give 

Pi(r)  =  cos(£(7r))  =  r*  +  C(r)(r  -  r*) 

where  now  C(r)  depends  on  the  bandwidth  parameter  WT.  The  Mean  Value 
Theorem  implies 

cos(7r)  -  cos(^(7r))  =  —  (7r  -  £(7r))sin(zr) 

where  zT  is  between  ~fr  and  E( jr).  If  we  define  yr  =  —  E(^T))  sin(zr), 

then  we  obtain  the  observed  contraction  mapping 

cos(7r)  =  r*  +  C(r)(r  -  r*)  +  yr  (22) 

This  expression  identifies  the  exact  error  incurred  by  replacing  cos(E(yr))  by 
cos(7r).  Now  consider  the  HK  algorithm 


rk+ i=cos(7rJ  (23) 

This  and  (  22)  give 

(24) 

rk+ 1  =  r*  +  [fic(  rj)](ri  r  )  d"  [  C{rff\yri  +  •  *  •  +  C{rk)yTk  i  yTk 

j= 1  3= 2 
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To  achieve  convergence,  the  bandwidth  must  go  to  0  at  a  certain  rate.  To 
uncover  this  rate,  we  resort  to  the  following  assumptions. 

Al.  Motivated  by  Proposition  3.1  and  Figures  1,2,  we  assume  that  in 
a  neighborhood  of  r*,  C{rk)  is  smaller  than  the  bandwidth  times  a  constant. 
Formally  we  assume  that  there  is  a  positive  constant  K  such  that 

0  <  C(rk)  <  KWk 

A2.  Motivated  by  Corollary  4.1,  we  assume  that 

\Vrk\  <  KWk 

where,  for  simplicity,  K  is  as  in  Al.. 

Remark  4.1.  Figures  1,2,  indicate  that  there  are  cases  where  C(rk)  is  very 
close  to  0  in  a  neighborhood  of  r*  in  addition  to  being  much  smaller  than 
the  bandwidth.  Thus,  Al  is  very  plausible.  Note  however,  that  we  stress 
the  requirement  of  being  in  a  neighborhood  of  r*,  for  otherwise,  when  the 
frequency  “is  not  captured,”  C(rk)  is  equal  to  1  (case  of  ideal  bandpass)  and 
convergence  toward  the  frequency  does  not  occur. 

Remark  4.2.  In  light  of  Corollary  4.1,  A2  is  very  plausible  also.  However, 
it  should  be  noted  that  in  Corollary  4.1,  the  parameter  r  is  fixed,  while  it  is 
random  in  A2. 


So,  under  these  assumptions, 

|rfc+i  -  r*|  <  KWk  {1  +  KWk-i  +  ■■■  +  ( KW2  •  •  •  KWk-i)} 


We  can  now  shrink  the  bandwidth  144,  as  k  — >  oo ,  at  a  rate  (for  example  at 
a  rate  of  1  /k)  that  guarantees  the  almost  sure  convergence  of  to  r*.  In 
this  case, 

7r*  -+  Wi 

almost  surely.  In  summary  we  have. 


n  KWj  (n  -  r*  +  1) 


(25) 


7=1 
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Theorem  4.1  Let  {Ztj  be  given  by  (  3),  and  let  {£r},  r  €  (  —  1,1),  be  a 
parametric  family  of  time  invariant  bandpass  linear  filters ,  with  bandwidth 
Wr,  for  which  the  fundamental  property  (  f)  holds,  and  one  that  passes  u>\ . 
Assume  that  the  filtered  process  {Zt(r)}  is  real  valued  for  all  r  €  (  —  1,1). 
Consider  the  recursion 

rk+1  =  cos(7rJ 

and  suppose  WTk  is  of  the  order  0(l/k).  If  Al  and  A2  hold,  then 


7 rk  ->  ^1 


(26) 


with  probability  one. 


4.1  Connection  between  C(r)  and  the  SNR 

The  signal  to  noise  ratio  (SNR)  for  the  process  (  3)  is  defined  by 

SNR  =  20  log10  —dB 


In  general,  the  SNR  is  defined  by, 


SNR  —  101og10 


Variance  of  the  filtered  signal 
Variance  of  the  filtered  noise 


dB 


from  which  it  follows  that  the  SNR  can  be  expressed  in  terms  of  the  contrac¬ 
tion  factor, 


SAMJ=101oglo{J--l}  (27) 

C{r) 

Thus,  the  SNR  and  C(r)  are  inversely  related: 

SNR ->oc  «  C(r)  — ►  0 
SNR  —*  — oo  C(r)  — >  1 

Because  we  can  control  yr  by  varying  the  bandwidth,  and  in  addition  in 
practice  yT  is  rather  small  (KS),  (  27)  shows  that  the  speed  of  convergence 
of  rk,  depends  primarily  on  the  contraction  factor  C(r).  The  smaller  C(r) 
is,  the  faster  is  the  speed  of  convergence. 
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4.2  Examples 

The  HK  algorithm  r*.+1  =  cos(7rjt)  was  applied  to  time  series  of  length 
N  =  20001  from  the  process  (  3)  with  Gaussian  white  noise  (t,  unless  stated 
otherwise,  using  the  real  counterpart  of  the  complex  filter  (  13), 


h(n;  r,  M )  = 


cos(n6(r)) / \/2M  +  1,  |n|  <  M 
0,  jnj  >  M 


(28) 


We  refer  to  this  filter  as  the  “cosine  filter”.  The  corresponding  squared  gain, 
\H(u;r,  M)\2,  is  closely  related  to  \H(uj-,r,  M)\2  in  (  14)  for  w  >  0, 

|7f  (W;  r,  A/)|2  =  i|  H(w,  r,  M)|2  +  0(~)  (29) 

Figure  3  shows  the  graph  of  \H(u;r,  M)\2,  u  >  0,  with  (Mi  =  2O,0(ri)  = 
1)  and  (M2  =  4O,0(r2)  =  1.1).  The  figure  illustrates  a  shift  in  location 
accompanied  by  a  reduced  bandwidth. 

As  can  be  seen  from  the  tables  in  the  following  examples,  the  convergence 
is  very  fast,  and  the  frequency  is  never  lost  when  Mk  increases  at  a  linear 
rate  such  as  {5, 15, 25, 35,  *  *  - }  -  Roughly,  this  is  equivalent  to  a  bandwidth 
that  shrinks  at  the  order  of  0(1/ k).  Also  note  that  the  sequence  of  observed 
zero-crossing  rates  {qv*}  is  monotone,  with  a  possible  exception  of  the  very 
beginning,  as  is  well  expected  from  (  7)  and  (  9). 

Very  similar  results  were  obtained  by  a  variety  of  other  bandpass  filters. 

It  should  be  noted  that  in  Y  and  KY  the  autocorrelation  was  estimated 
from  the  sample  autocorrelation  of  the  filtered  processes,  while  here  the  es¬ 
timates  are  obtained  from  the  observed  zero-crossing  rate  after  filtering  as 
suggested  from  HK  and  KS.  Evidently,  the  two  methods  of  estimation,  in 
conjunction  with  the  CM  method,  give  very  similar  frequency  estimates  for 
sufficiently  long  data  sets. 

lMore  precisely,  for  a  given  M*,  the  zero-crossing  count  is  computed  from  time  series 
of  length  2000  +  2 Mk- 
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4.2.1  Example  1:  A  simple  case 

ui  =  0.8,  0(ro)  =  0.4,  Ai  =  —0.91333,  Bi  =  0.06975,  SNR  =  0 dB, 


k 

Mk 

rk 

lrk 

0 

5 

0.92106 

0.65221 

1 

15 

0.79475 

0.75279 

2 

25 

0.72979 

0.79994 

3 

35 

0.69675 

0.79994 

4 

45 

0.69675 

0.79994 

5 

55 

0.69675 

0.79994 

Cj\  =  0.79994 

For  comparison,  the  five  highest  periodogram  ordinates,  ordered  by  mag¬ 
nitude,  were  at  the  frequencies: 

0.80185,  0.80485,  2.26793,  0.74201,  0.83776 
with  the  maximum  occurring  at  0.80185. 

4.2.2  Example  2:  NonGaussian  and  colored  noise 

When  the  fundamental  property  (  4)  holds  for  white  noise,  it  also  holds  ap¬ 
proximately  for  slowly  varying  continuous  spectral  densities,  as  long  as  the 
bandwidth  is  sufficiently  narrow.  In  this  respect,  it  is  of  interest  to  apply 
the  HK  algorithm  with  the  cosine  filter  and  colored  Gaussian  noise.  It  is 
also  interesting  to  run  the  algorithm  when  the  noise  is  white  but  nonGaus- 
sian.  Accordingly,  we  have  applied  the  algorithm  when  (t  is  Gaussian  white 
noise  (GWN),  lognormal  white  noise  (LNWN)  (with  parameters  0,1  ,  and 
properly  centered  and  scaled),  and  a  first-order  moving  average  (MA(1)), 
ut  =  Ct  ~  O.8G-1,  and  Gaussian  white  noise.  In  all  three  cases  the  si¬ 
nusoidal  component  was  identical.  The  results  are  given  in  the  following 
table,  and  show  identical  convergence.  Evidently,  in  the  bandpass  case,  the 
algorithm  is  quite  robust. 
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Ui  =  0.71,  6(r0)  =  1.1,  A\  -  —0.954817,  B,  =  1.00687, 
SNR  =  0  dB, 


GWN _ LNWN _  MA(1) 


B 

Mk 

Irk 

lrk 

7r* 

0 

5 

0.993240 

0.913089 

0.997955 

i 

15 

0.920947 

0.897373 

0.928805 

2 

25 

0.812508 

0.766932 

0.848654 

3 

35 

0.766932 

0.710355 

0.748073 

4 

45 

0.719785 

0.710355 

0.710355 

5 

55 

0.710355 

0.710355 

0.710355 

6 

65 

0.710355 

0.710355 

0.710355 

loi  =  0.710355 

ib,  =  0.710355 

=  0.710355 
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4.2.3  Example  3:  Mixed  spectrum  noise 

In  this  example  we  compare  the  convergence  of  the  HK  algorithm  in  two 
cases  when  the  noise  (t  is  Gaussian  white  noise,  and  then  when  it  is  Gaussian 
white  noise  plus  a  random  Gaussian  sinusoid  (mixed  spectrum  noise)  with 
frequency  2.5.  In  both  cases  the  “signal”  is  a  random  Gaussian  sinusoid  with 
frequency  uq  =  1.25,  and  the  search  starts  at  0(ro)  =  0.5.  We  have:  A1  == 
0.30041,  Bx  =  —1.44635,  the  amplitude  of  the  sinusoidal  noise  component  is 
0.947085,  and  SNR  =  OdB.  Again,  from  the  table,  the  convergence  is  very 
satisfactory. 

This  example  shows  that  our  algorithm  can  be  applied  separately  to  each 
frequency  using  bandpass  filters.  In  fact  when  #(r0)  =  2.3  the  7rjfc  converges 
to  2.50039,  and  the  other  frequency  is  detected  as  expected. 


GWN  Mixed  Spectrum  Noise 


k 

Mk 

7r* 

7r* 

0 

5 

0.93509 

1.09854 

1 

15 

1.05296 

1.23684 

2 

25 

1.16140 

1.24941 

3 

35 

1.16454 

1.24941 

4 

45 

1.22426 

1.24941 

5 

55 

1.24941 

1.24941 

6 

65 

1.24941 

1.24941 

7 

75 

1.24941 

1.24941 

8 

85 

1.24941 

1.24941 

ux  =  1.24941 

«!  =  1.24941 

With  Gaussian  white  noise,  the  highest  peaks  in  the  periodogram,  ordered 
by  magnitude,  occur  at  the  frequencies, 

1.25366, 1.25068, 1.18815, 1.09286, 1.30726 

with  the  highest  peak  occurring  at  1.25366.  This  is  the  periodogram  estimate. 
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4.2.4  Example  4:  Moderately  low  SNR 

u\  =  0.4,  0(ro )  =  0.2,  Ax  =  —0.537689,  B\  =  -0.423486,  SNR  =  -6.0206dR, 


k 

Mk 

rk 

7r* 

0 

5 

0.98007 

0.33160 

1 

15 

0.94552 

0.40075 

2 

25 

0.92077 

0.40075 

3 

35 

0.92077 

0.40075 

4 

45 

0.92077 

0.40075 

5 

55 

0.92077 

. 

0.40075 

W!  =  0.40075 

The  highest  periodogram  ordinates  correspond  to  the  frequencies 
0.40392, 1.93283, 0.46376, 2.01660, 2.05550 
The  largest  peak  occurs  at  0.40392. 
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4.2.5  Example  5:  low  SNR 

ui  =  1.95,  0(ro)  =  1.8,  =  0.174937,  Bx  =  -0.263656,  SNR  =  -13.9794rfR. 


k 

Mk 

Tk 

lrk 

0 

5 

-0.22720 

1.87647 

1 

15 

-0.30094 

1.89690 

2 

25 

-0.32035 

1.91576 

3 

35 

-0.33816 

1.93776 

4 

45 

-0.35878 

1.95191 

5 

55 

-0.37195 

1.94876 

6 

65 

-0.36903 

1.94876 

7 

75 

-0.36903 

1.94876 

8 

85 

-0.36903 

1.94876 

ui  =  1.94876 

The  five  highest  peaks  in  the  periodogram  are  at 

2.18584, 2.36873, 1.27729, 0.51327, 1.02655 
where  the  highest  power  is  given  to  2.18584. 

4.2.6  Example  6:  Detection  of  a  diurnal  cycle 

GATE  stands  for  GARP  (Global  Atmospheric  Research  Program)  Atlantic 
Tropical  Experiment  (see  Simpson  (1988)).  The  experiment,  conducted  in 
1974  in  the  eastern  Atlantic  off  the  coast  of  west  Africa,  produced  large 
snapshots  of  instantaneous  tropical  rain  rate  collected  in  15  minutes  intervals. 
The  experiment  was  conducted  in  stages  of  which  the  first  is  referred  to  as 
GATE  I.  The  hourly  GATE  I  data2,  averaged  over  a  region  of  280  x  280  km2, 
are  plotted  in  Figure  4. 

2Thanks  are  due  to  Dr.  Tom  Bell,  NASA,  Goddard  SFC,  for  kindly  supplying  the 
hourly  GATE  data. 
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From  an  atmospheric  dynamics  point  of  view,  it  is  useful  to  know  whether 
tropical  rainfall  contains  a  significant  diurnal  cycle  (Simpson  (1988),  p.  40). 
It  is  believed  that  the  GATE  I  data  set  does  contain  a  significant  diurnal 
cycle.  However,  so  far,  the  detection  of  the  diurnal  cycle  in  GATE  has  not 
been  quite  decisive. 

We  have  applied  the  HK  algorithm,  using  the  cosine  filter,  to  the  centered 
hourly  GATE  I  data  in  search  of  a  diurnal  cycle.  As  can  be  seen  from  the 
following  table,  starting  from  various  points,  the  convergence  is  toward  the 
value  0.262675  radians  per  hour.  This  gives  a  period  of 


2tt 

0.262675 


=  23.92  hours 


which  gives  more  credence  to  the  the  existence  of  a  diurnal  cycle  in  GATE  I. 


fl(r0)  =  0.50  fl(r0)  =  0.45  fl(r0)  =  0.30 


k 

Mk 

7r* 

7r* 

7  Tk 

0 

55 

0.483322 

0.441294 

0.273182 

1 

56 

0.420280 

0.441294 

0.262675 

2 

57 

0.367745 

0.399266 

0.262675 

3 

58 

0.346731 

0.346731 

0.262675 

4 

59 

0.304703 

0.304703 

0.262675 

5 

60 

0.294196 

0.294196 

0.262675 

6 

61 

0.283689 

0.283689 

0.262675 

7 

62 

0.262675 

0.262675 

0.262675 

8 

63 

0.262675 

0.262675 

0.262675 

9 

64 

0.262675 

0.262675 

* 

0.262675 

wt  =  0.262675 

Cjx  =  0.262675 

wi  =  0.262675 
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4.2.7  Example  7:  Comparison  between  two  filters 

The  real  counterpart  of  the  complex  filter  defined  by  (  19)  is  obtained  from 
the  real  impulse  response, 


h(n ;  r,  M)  =  < 


M 

n 

0, 


cos (0(r)rc),  n  =  0, 

otherwise 


(30) 


We  refer  to  this  filter  as  the  “cosine-binomial”  filter.  To  achieve  a  degree  of 
precision  comparable  to  the  one  achieved  by  the  cosine  filter  (  28),  usually  a 
larger  M  is  needed  for  the  cosine-binomial  filter  (  30). 

The  following  table  compares  the  performance  of  the  cosine  filter  with 
that  of  the  cosine-binomial  filter  considering  only  six  iterations  of  the  HK 
algorithm.  Here  ux  =  1.0,  0(ro )  =  0.6,  Ax  =  — 0.181463,  -Si  =  —1.01382, 
SNR  =  0  dB. 


Cosine  filter  Cos-bin  filter 


k 

Mk 

7rk 

7r* 

0 

5 

0.89580 

1.04667 

1 

10 

0.99953 

1.04039 

2 

15 

0.99953 

1.03096 

3 

20 

1.00110 

1.01210 

4 

25 

1.00110 

1.00581 

5 

30 

1.00110 

1.00424 

ux  =  1.00110 

Cj\  -  1.00424 
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5  Summary 

Our  observation  is  that  for  sufficiently  long  data  records  (e.g.  N  >  1000), 
and  for  a  reasonable  signal  to  noise  ratio  (e.g.  SNR  >  0),  the  original 
CM  algorithm,  as  given  in  HK  in  terms  of  HOC  from  zero-crossings,  per¬ 
forms  remarkably  well  when  the  method  employs  bandpass  filters.  Further¬ 
more,  the  method  requires  O(N)  computational  complexity,  compared  with 
0(N  log  N)  complexity  associated  with  FFT  based  methods.  In  almost  all 
the  cases  considered  in  the  above  examples,  the  HK  algorithm  outperformed 
the  periodogram.  This  has  been  the  case  with  numerous  other  examples, 
using  a  variety  of  filters.  However,  combining  both  methods  may  result  in 
the  improvement  of  both. 

We  emphasized  the  single  frequency  case.  However,  an  extension  to  the 
multiple  frequency  case  is  rather  straightforward  because  the  CM  method, 
when  employing  bandpass  filters,  can  be  applied  independently  to  nonover¬ 
lapping  frequency  bands  as  done  in  KY.  That  this  is  indeed  feasible  was 
already  seen  in  Example  3.  As  such,  this  scheme  is  ideal  for  parallel  process¬ 
ing. 

We  have  used  only  a  few  filters  in  connection  with  the  CM  method,  but 
of  course  many  more  are  available.  Results  concerning  different  filters  will 
be  dealt  with  elsewhere. 

There  is  one  more  thing.  We  have  mentioned  earlier  that  for  a  zero- 
mean  stationary  Gaussian  process  with  a  mixed  spectrum,  the  observed  zero- 
crossing  rate  is  not  necessarily  consistent  for  its  expectation.  More  precisely, 
in  KS  it  is  shown  that  if  the  spectrum  contains  two  or  more  atoms  (that 
is,  two  or  more  jumps  in  the  spectral  distribution  function),  then  qy  does 
not  converge  to  a  constant  as  A'’  — >  oo.  However,  when  the  spectrum  is 
a  mixture  of  an  atomic  component  with  a  single  atom,  and  an  absolutely 
continuous  component  (noise),  %  is  almost  surely  consistent  provided  both 
the  sinusoidal  component  and  the  noise  component  have  the  exact  same  first 
order  autocorrelation.  This,  rather  esoteric,  case  is  precisely  what  the  HK 
algorithm  gives  in  the  limit.  To  see  this,  consider  the  process  (  3),  and  let 

Xt  =  A\  cos(u;i t)  -f-  B\  sin(u;if) 


Then 


zt  =  xt  +  Ct 
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Let  {Cr}  satisfy  (  4).  Then  for  r  =  r*  =  cos(u7),  we  have, 


cAZ)t  =  cAX)t  +  cT.{Qt 

Because  r*  is  a  fixed  point  of  Pi(r),  £r.( X)t  and  £r.(£)<  have  r*  as  their  first 
order  autocorrelation,  and  so 


7r*  E( 7r.)  =  u>x 

almost  surely.  This  means  that  when  the  HK  algorithm  converges  (bandpass 
filters  or  not),  the  limiting  zero-crossing  rate  is  almost  surly  a  constant. 

Acknowledgement.  The  authors  are  grateful  to  Eric  Slud  for  helpful  com¬ 
ments. 

Figures 


Figure  1.  (a).  An  attracting  fixed  point  that  practically  coincides  with 
cos(0.8)  in  pj(r,  20),  ui\  —  0.8,  M  =  20,  from  the  complex  exponential  filter. 
The  derivative  of  px(a,  20)  in  a  neighborhood  of  cos(0.8)  is  very  close  to  0. 
(b).  The  squared  gain  (  14)  with  M  =  20,  and  centered  at  0  =  0.8. 

Figure  2.  (a).  An  attracting  fixed  point  that  practically  coincides  with 
cos(0.8)  in  pi(a,  100),  =  0.8,  M  =  100,  from  the  complex  filter  in  HK. 

The  derivative  of  pi(a,  100)  in  a  neighborhood  of  cos(0.8)  is  close  to  0.  (b)„ 
The  normalized  squared  gain  (  20)  with  M  =  100,  and  centered  at  0  —  0.8. 

Figure  3.  The  squared  gain  of  the  cosine  filter  centered  at  6  —  1,  with 
M  =  20,  and  centered  at  9  =  1.1  with  M  =  40. 

Figure  4.  GATE  I  data:  interpolated  hourly  rain  rate  averaged  over  a 
280  x  280  km 2  area  in  the  eastern  Atlantic.  N  =  450. 
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Figure  1.  (a).  An  attracting  fixed  point  that  practically  coincides  with 
cos(0.8)  in  pi(r,  20),  ui\  ~  0.8,  M  =  20,  from  the  complex  exponential  filter. 
The  derivative  of  />!(a,20)  in  a  neighborhood  of  cos(0.8)  is  very  close  to  0. 
(b).  The  squared  gain  (  14)  with  M  =  20,  and  centerd  at  0  =  0.8. 


ro(r)  From  Fejer  Kernel:  M=20,  cos (0 . 8) =0 . 697 
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Figure  2.  (a).  An  attracting  fixed  point  that  practically  coincides  with 
cos(0.8)  in  pi(a,  100),  uq  =  0.8,  M  =  100,  from  the  complex  filter  in  HI\. 
The  derivative  of  pi(a,  100)  in  a  neighborhood  of  cos(0.8)  is  close  to  0.  (b). 
The  normalized  squared  gain  (  20)  with  M  =  100,  and  centerd  at  B  =  0.8. 


ro(r)  From  Complex  Filter:  M=100,  cos (0 . 8) =0. 697 


[cos (w-Theta) /2]~2M,  M=100,  Theta=0.8 


Figure  3.  The  squared  gain  of  the  cosine  filter  centered  at  6 
M  —  20,  and  centered  at  9  =  1.1  with  M  =  40. 
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Cosine  Kernel:  (M=20, w ( 1) =1 ) ,  (M=40, w ( 1 ) =1 . 1 ) 


Figure  4.  GATE  I  data:  interpolated  hourly  rain  rate  averaged 
280  x  280  km 2  area  in  the  eastern  Atlantic.  N  =  450. 
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